"""
特征提取与降维案例：PCA、LDA、t-SNE
演示各种降维方法及其应用
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
from sklearn.datasets import load_iris, load_digits
from sklearn.decomposition import PCA, KernelPCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
import warnings
warnings.filterwarnings('ignore')

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 80)
print("特征提取与降维案例：PCA、LDA、t-SNE")
print("=" * 80)

# 1. 加载数据
print("\n1. 加载数据集")
print("=" * 80)

# 使用鸢尾花数据集
iris = load_iris()
X_iris = iris.data
y_iris = iris.target
feature_names = iris.feature_names
target_names = iris.target_names

print(f"鸢尾花数据集形状：{X_iris.shape}")
print(f"特征数量：{X_iris.shape[1]}")
print(f"样本数量：{X_iris.shape[0]}")
print(f"类别：{target_names}")
print(f"特征名称：{list(feature_names)}")

# 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_iris)

# ============================================================================
# PCA - 主成分分析
# ============================================================================

print("\n" + "=" * 80)
print("2. PCA - 主成分分析 (Principal Component Analysis)")
print("=" * 80)
print("原理：找到数据方差最大的方向，进行线性变换")
print("特点：无监督、线性降维")
print("目标：最大化投影后的方差")

# 2.1 应用PCA
pca = PCA()
X_pca = pca.fit_transform(X_scaled)

# 查看解释方差比
explained_variance_ratio = pca.explained_variance_ratio_
cumulative_variance_ratio = np.cumsum(explained_variance_ratio)

print("\n各主成分的解释方差比：")
for i, (var, cum_var) in enumerate(zip(explained_variance_ratio, cumulative_variance_ratio)):
    print(f"PC{i+1}: {var:.4f} (累计: {cum_var:.4f})")

# 2.2 选择主成分数量
print("\n选择主成分数量的策略：")
print("1. 累计解释方差达到85%-95%")
print("2. 碎石图（Scree Plot）找拐点")
print("3. Kaiser准则：特征值>1")

# 使用2个主成分
pca_2d = PCA(n_components=2)
X_pca_2d = pca_2d.fit_transform(X_scaled)

print(f"\n降维到2维后的解释方差比: {pca_2d.explained_variance_ratio_}")
print(f"累计解释方差: {sum(pca_2d.explained_variance_ratio_):.4f}")

# 2.3 主成分的含义
print("\n主成分的组成（载荷矩阵）：")
components_df = pd.DataFrame(
    pca_2d.components_,
    columns=feature_names,
    index=['PC1', 'PC2']
)
print(components_df)

# 2.4 PCA可视化
fig = plt.figure(figsize=(18, 5))

# 碎石图
ax1 = fig.add_subplot(131)
ax1.plot(range(1, len(explained_variance_ratio) + 1), 
         explained_variance_ratio, 'bo-', label='单个方差')
ax1.plot(range(1, len(cumulative_variance_ratio) + 1), 
         cumulative_variance_ratio, 'rs-', label='累计方差')
ax1.axhline(y=0.95, color='g', linestyle='--', label='95%阈值')
ax1.set_xlabel('主成分')
ax1.set_ylabel('解释方差比')
ax1.set_title('PCA碎石图')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2D散点图
ax2 = fig.add_subplot(132)
colors = ['red', 'green', 'blue']
for i, target_name in enumerate(target_names):
    mask = y_iris == i
    ax2.scatter(X_pca_2d[mask, 0], X_pca_2d[mask, 1], 
                c=colors[i], label=target_name, alpha=0.6, s=50)
ax2.set_xlabel(f'PC1 ({pca_2d.explained_variance_ratio_[0]:.2%})')
ax2.set_ylabel(f'PC2 ({pca_2d.explained_variance_ratio_[1]:.2%})')
ax2.set_title('PCA降维到2D')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 载荷图
ax3 = fig.add_subplot(133)
for i, feature in enumerate(feature_names):
    ax3.arrow(0, 0, components_df.loc['PC1', feature], 
              components_df.loc['PC2', feature],
              head_width=0.05, head_length=0.05, fc='blue', ec='blue')
    ax3.text(components_df.loc['PC1', feature] * 1.15,
             components_df.loc['PC2', feature] * 1.15,
             feature, fontsize=9)
ax3.set_xlabel('PC1')
ax3.set_ylabel('PC2')
ax3.set_title('PCA载荷图')
ax3.grid(True, alpha=0.3)
ax3.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax3.axvline(x=0, color='k', linestyle='-', linewidth=0.5)

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/5-特征提取与降维/PCA分析.png',
            dpi=300, bbox_inches='tight')
print("\nPCA可视化已保存")

# ============================================================================
# LDA - 线性判别分析
# ============================================================================

print("\n" + "=" * 80)
print("3. LDA - 线性判别分析 (Linear Discriminant Analysis)")
print("=" * 80)
print("原理：最大化类间距离，最小化类内距离")
print("特点：有监督、线性降维")
print("目标：最大化类别可分性")

# 应用LDA
lda = LDA(n_components=2)
X_lda = lda.fit_transform(X_scaled, y_iris)

print(f"\nLDA降维到2维")
print(f"解释方差比: {lda.explained_variance_ratio_}")
print(f"累计解释方差: {sum(lda.explained_variance_ratio_):.4f}")

# LDA可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# LDA 2D散点图
for i, target_name in enumerate(target_names):
    mask = y_iris == i
    axes[0].scatter(X_lda[mask, 0], X_lda[mask, 1],
                    c=colors[i], label=target_name, alpha=0.6, s=50)
axes[0].set_xlabel(f'LD1 ({lda.explained_variance_ratio_[0]:.2%})')
axes[0].set_ylabel(f'LD2 ({lda.explained_variance_ratio_[1]:.2%})')
axes[0].set_title('LDA降维到2D')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# PCA vs LDA对比
axes[1].scatter(X_pca_2d[:, 0], X_pca_2d[:, 1], 
                c=y_iris, cmap='viridis', alpha=0.3, s=30, label='PCA')
axes[1].scatter(X_lda[:, 0], X_lda[:, 1],
                c=y_iris, cmap='plasma', alpha=0.6, s=30, marker='^', label='LDA')
axes[1].set_xlabel('第一维度')
axes[1].set_ylabel('第二维度')
axes[1].set_title('PCA vs LDA对比')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/5-特征提取与降维/LDA分析.png',
            dpi=300, bbox_inches='tight')
print("LDA可视化已保存")

# ============================================================================
# t-SNE - t分布随机邻域嵌入
# ============================================================================

print("\n" + "=" * 80)
print("4. t-SNE - t分布随机邻域嵌入")
print("=" * 80)
print("原理：保持数据的局部结构，非线性降维")
print("特点：无监督、非线性、主要用于可视化")
print("注意：计算成本高，结果受参数影响大")

# 使用手写数字数据集（更适合t-SNE）
digits = load_digits()
X_digits = digits.data
y_digits = digits.target

print(f"\n手写数字数据集形状：{X_digits.shape}")
print(f"特征数量：{X_digits.shape[1]}")

# 标准化
X_digits_scaled = StandardScaler().fit_transform(X_digits)

# 应用t-SNE（使用部分数据以加快速度）
n_samples = 1000
indices = np.random.choice(len(X_digits_scaled), n_samples, replace=False)
X_subset = X_digits_scaled[indices]
y_subset = y_digits[indices]

print(f"\n使用{n_samples}个样本进行t-SNE...")

# 不同困惑度的t-SNE
perplexities = [5, 30, 50]
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for idx, perplexity in enumerate(perplexities):
    print(f"困惑度 = {perplexity}")
    tsne = TSNE(n_components=2, perplexity=perplexity, random_state=42)
    X_tsne = tsne.fit_transform(X_subset)
    
    scatter = axes[idx].scatter(X_tsne[:, 0], X_tsne[:, 1],
                                c=y_subset, cmap='tab10', alpha=0.6, s=20)
    axes[idx].set_title(f't-SNE (困惑度={perplexity})')
    axes[idx].set_xlabel('t-SNE 1')
    axes[idx].set_ylabel('t-SNE 2')
    plt.colorbar(scatter, ax=axes[idx])

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/5-特征提取与降维/tSNE分析.png',
            dpi=300, bbox_inches='tight')
print("t-SNE可视化已保存")

# ============================================================================
# 核PCA - 非线性PCA
# ============================================================================

print("\n" + "=" * 80)
print("5. 核PCA (Kernel PCA)")
print("=" * 80)
print("原理：使用核技巧进行非线性降维")
print("特点：可以捕捉非线性关系")

# 应用不同核函数的PCA
kernels = ['linear', 'rbf', 'poly']
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for idx, kernel in enumerate(kernels):
    kpca = KernelPCA(n_components=2, kernel=kernel, random_state=42)
    X_kpca = kpca.fit_transform(X_scaled)
    
    for i, target_name in enumerate(target_names):
        mask = y_iris == i
        axes[idx].scatter(X_kpca[mask, 0], X_kpca[mask, 1],
                         c=colors[i], label=target_name, alpha=0.6, s=50)
    axes[idx].set_xlabel('KPC1')
    axes[idx].set_ylabel('KPC2')
    axes[idx].set_title(f'核PCA ({kernel}核)')
    axes[idx].legend()
    axes[idx].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/5-特征提取与降维/核PCA分析.png',
            dpi=300, bbox_inches='tight')
print("核PCA可视化已保存")

# ============================================================================
# 降维对模型性能的影响
# ============================================================================

print("\n" + "=" * 80)
print("6. 降维对模型性能的影响")
print("=" * 80)

# 分割数据
X_train, X_test, y_train, y_test = train_test_split(
    X_scaled, y_iris, test_size=0.3, random_state=42
)

results = {}

# 原始数据
rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)
results['原始数据(4维)'] = accuracy_score(y_test, rf.predict(X_test))

# 不同维度的PCA
for n_components in [2, 3]:
    pca_temp = PCA(n_components=n_components)
    X_train_pca = pca_temp.fit_transform(X_train)
    X_test_pca = pca_temp.transform(X_test)
    
    rf = RandomForestClassifier(n_estimators=100, random_state=42)
    rf.fit(X_train_pca, y_train)
    results[f'PCA({n_components}维)'] = accuracy_score(y_test, rf.predict(X_test_pca))

# LDA
lda_temp = LDA(n_components=2)
X_train_lda = lda_temp.fit_transform(X_train, y_train)
X_test_lda = lda_temp.transform(X_test)

rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X_train_lda, y_train)
results['LDA(2维)'] = accuracy_score(y_test, rf.predict(X_test_lda))

print("\n不同降维方法的模型准确率：")
for method, acc in results.items():
    print(f"{method}: {acc:.4f}")

# 可视化结果
fig, ax = plt.subplots(figsize=(10, 6))
methods = list(results.keys())
accuracies = list(results.values())
colors_bar = ['gray', 'skyblue', 'lightblue', 'lightgreen']

bars = ax.barh(methods, accuracies, color=colors_bar)
ax.set_xlabel('准确率')
ax.set_title('不同降维方法对模型性能的影响')
ax.set_xlim([0.8, 1.0])
ax.grid(True, alpha=0.3, axis='x')

for i, (bar, acc) in enumerate(zip(bars, accuracies)):
    ax.text(acc, i, f' {acc:.4f}', va='center')

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/5-特征提取与降维/降维性能对比.png',
            dpi=300, bbox_inches='tight')
print("\n性能对比可视化已保存")

# ============================================================================
# 降维方法对比总结
# ============================================================================

print("\n" + "=" * 80)
print("7. 降维方法对比总结")
print("=" * 80)
print("""
┌─────────────┬──────────────┬──────────────┬─────────────────┬──────────────────┐
│   方法      │   类型       │   是否监督   │    适用场景     │      特点        │
├─────────────┼──────────────┼──────────────┼─────────────────┼──────────────────┤
│ PCA         │ 线性         │ 无监督       │ 去相关、降噪    │ 最大化方差       │
│             │              │              │ 数据压缩        │ 快速、可解释     │
├─────────────┼──────────────┼──────────────┼─────────────────┼──────────────────┤
│ LDA         │ 线性         │ 有监督       │ 分类任务        │ 最大化可分性     │
│             │              │              │ 特征提取        │ 类别数-1维上限   │
├─────────────┼──────────────┼──────────────┼─────────────────┼──────────────────┤
│ 核PCA       │ 非线性       │ 无监督       │ 非线性关系      │ 灵活、计算复杂   │
│ Kernel PCA  │              │              │                 │                  │
├─────────────┼──────────────┼──────────────┼─────────────────┼──────────────────┤
│ t-SNE       │ 非线性       │ 无监督       │ 数据可视化      │ 保持局部结构     │
│             │              │              │ 聚类分析        │ 计算慢、难调参   │
└─────────────┴──────────────┴──────────────┴─────────────────┴──────────────────┘

选择建议：

1. **PCA**
   - 数据预处理、去除相关性
   - 需要可解释的主成分
   - 数据近似线性关系
   - 快速降维

2. **LDA**
   - 分类任务的特征提取
   - 需要最大化类别区分度
   - 有标签数据
   - 类别数较少

3. **核PCA**
   - 数据有明显非线性关系
   - PCA效果不好时
   - 计算资源充足

4. **t-SNE**
   - 高维数据可视化（降到2D/3D）
   - 探索数据聚类结构
   - 不用于特征工程（仅可视化）
   - 结果不用于后续建模

参数选择：

PCA:
- n_components: 累计解释方差85%-95%，或碎石图拐点

LDA:
- n_components: 最多为类别数-1

t-SNE:
- perplexity: 5-50，数据量大用大值
- learning_rate: 10-1000
- n_iter: 至少1000

注意事项：
1. 降维前要标准化数据
2. 在训练集上fit，在测试集上transform
3. 保存降维模型用于新数据
4. 降维可能损失信息，需权衡
5. t-SNE每次运行结果不同，不可复现用于预测
""")

# 8. 实际应用示例
print("\n" + "=" * 80)
print("8. 实际应用Pipeline示例")
print("=" * 80)

from sklearn.pipeline import Pipeline

# 创建包含降维的Pipeline
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('pca', PCA(n_components=2)),
    ('classifier', RandomForestClassifier(n_estimators=100, random_state=42))
])

# 训练
X_train_raw, X_test_raw, y_train, y_test = train_test_split(
    X_iris, y_iris, test_size=0.3, random_state=42
)

pipeline.fit(X_train_raw, y_train)
score = pipeline.score(X_test_raw, y_test)

print(f"\nPipeline准确率: {score:.4f}")
print("\nPipeline步骤：")
print("1. StandardScaler: 标准化")
print("2. PCA: 降维到2维")
print("3. RandomForest: 分类")

# 保存Pipeline
import joblib
joblib.dump(pipeline, '/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/5-特征提取与降维/pipeline.pkl')
print("\nPipeline已保存")

print("\n" + "=" * 80)
print("案例完成！")
print("=" * 80)
